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Abstract 

Monthly averages of the sunspot number visible on the sun, observed from 1749, Zurich 
Observatory, and from 1848 other observatories, have been analyzed. This time signal 
presents a frequency power spectra with a clear 1//" behavior with a ~ 0.8±0.2. The well 
known cycle of approximately 11 years, clearly present in the spectrum, does not produce 
a sensible distortion of that behavior. The eventual characterization of the sunspot time 
series as a fractal is analyzed by means of the detrended fluctuation analysis. The jump 
size distribution of the signal is also studied. 



1 Introduction 

A sunspot [1] is a relatively cool area on the solar face that appears as a dark blemish. The 
number of sunspots is continuously changing in time in a random fashion and constitutes 
a typical, a priori, random time series. Because of the symmetry of the twisted magnetic 
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lines that originate sunspots, they are generally seen in pairs or in groups of pairs at both 
sides of the solar equator. As the sunspot cycle progresses, spots appear closer to the 
sun's equator giving rise to the so called "butterfly diagram" in the latitude distribution. 
Magnetic fields about the sunspots are very strong and keep heat out of these regions 
on the sun surface. They are formed when magnetic field lines are twisted and poke 
through the solar photosphere. The twisted magnetic fields above sunspots are sites where 
solar flares are observed. It has been found that chromospheric flares show a very close 
statistical relationship with sunspots. This is illustrated by the Waldmeir relationship 
[1] between the mean number of flares per day, E, and the sunspot number, R, namely 
E — 0.061 R. Being the number of sunspots much larger than the corresponding to flares, 
it seems more appropriate for statistical analysis. Moreover, this connection between 
sunspots and solar magnetic field is particularly motivating for a further analysis of their 
frequency of appearance. In summary, being a consequence of the solar cycle, the variation 
of the number of sunspots is closely related to the so-called solar activity. In fact, the 
solar maximum activity corresponds to many sunspots present, while the minimum of 
activity is related to very few sunspots. 

Another effect of the solar activity is the presence of fluctuations in the interplanetary 
magnetic fleld. Experimental data suggest [2] a 1// dependence for the frequency spectra 
of its fluctuations in the range of large times, namely, 2.7 x 10~^ Hz to 1.8 x 10^^ Hz. In 
this analysis it is argued that on the basis of the observed behavior there is a superpo- 
sition of signals with scale invariant distributions of correlation times. In other words, a 
kinematical superposition of signals due to scale- invariant reconnection of magnetic struc- 
tures, near the solar surface, gives rise to the 1/f spectrum. It should be also noticed that 
at higher frequencies the interplanetary spectrum is known to have approximatively an 
behavior; this time associated with magnetohydrodynamics (MHD) homogeneous 
turbulence. 

The result of our analysis of the sunspot number distribution along almost 250 years 
shows a clear correlation among all these phenomena. 

After introducing the data and the different methods of analysis used, we present the 
results with a further comment on their implications. 

2 Sunspot Data 

There is experimental evidence and data registers on sunspots since Gahleo (1610). Af- 
terwards, the sunspot data were provided by the Ziirich Observatory from 1749 to 1981. 
Since January 1981, the data are provided by the Sunspot Index Data Center (SIDC) [3]. 

In Figure 1 the monthly measured number of sunspots is plotted versus time. This 
plot includes 44,444 data points. The approximate 11 years cycle can clearly be seen 
from this plot. 
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Figure 1: Observed spots number as a function of time 

3 Analysis of the Time Series 

The analysis of a given temporal signal with apparently random fluctuations starts by 
asking if the value of the signal in a given instant has any correlation with the signal 
in a later time. The standard statistical tool for describing the signal is the temporal 
correlation function and the corresponding power spectrum or frequency spectrum also 
called spectral density. In general, random time dependent perturbations give rise to 
random noise characterized by a frequency spectrum following a power law 1//^ with 
—2 < (3 < 0. Such a power law is found in a rich variety of physical systems exhibiting 
some kind of episodic activity. When (3 —>■ 1, the correlation times grow strongly. This 
behavior is closely related to the concept of fractal and self-similarity. In many cases, the 
dynamics of the corresponding physical system is attracted towards a critical state and 
the system is said to present self-organized criticality. 

The concepts of fractal and self- similarity are directly related [4]. Fractals are self- 
similar under an isotropic reduction or scale transformation. 

In many cases of interest, the concept of self-affine fractality or self-affinity has shown 
to be very fruitful. In this "weak" form of similarity, the sequence of data looks alike 
when the reduction of the time scale is made with a different scaling factor than the one 
used for the measured variable. This means invariance under anisotropic scale changes. 
The random walk or Brownian movement defined by a function B{t) that determines 
the instant position of the particle is a well-known example. If -B(O) = 0, then at any 
time t the characteristic function of the random walk verifies B{fit) for any fi 0. 

The Brownian exponent is precisely Hi = 1/2. In any case it is possible to characterize 
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the fractal behavior of a system through a non-integer dimension usually called fractal 
dimension related to the above mentioned Hi. 

An alternative characterization of fractal time series is by means of the roughness or 
Hurst exponent Hu = 1 + Hi that measures the persistency of the fluctuations related to 
the time series [5] . 

There have been different proposals for measuring these exponents. One of them is 
the analysis of the variance of increments at a given lag, also known with the geophysical 
name of semivariogram [6]. The detrcndcd fluctuation analysis DFA [8] has shown to be 
an alternative procedure, some times a better one, for analyzing fluctuations in time series 
and a way of defining another similar exponent Ha- There are other possibihties as the 
mobile average analysis introduced in reference [9]. 

The statistical properties of the sunspot number fluctuations, were analyzed by means 
of the different tools just mentioned, plus the jump size distribution. 



3.1 Analysis of the Frequency Power Spectrum 

The analysis of correlation in time series starts with the temporal correlation function 
deflned by 

G(r) = (x(to) x{to + r))t„ - {x{to))l (1) 

Clearly, whenever no temporal correlation is present, one has G{t) = 0. On the other 
hand, the rate with which G{t) decreases from the average value of the instantaneous 
fluctuation towards zero measures the correlation time that evaluates the memory effects 
present in the signal. 

The frequency spectrum is deflned in terms of the amplitude of the Fourier transform 
of the time signal squared, namely 

2 

(2) 

In the particular case of a stationary process, i.e. its behavior is independent of the 
particular instant chosen to start the observation, the last expression ends in the cosine 
Fourier transform 

S{f)^2 dr G{t) cos{2TTf t) (3) 

^0 

In general, a random time dependent fluctuation gives rise to a random noise characterized 
by a frequency spectrum following a power law 

S{f) oc (4) 

over a large range of frequencies, with well known examples in the range —2 < /3 < 0. 

There exists a heuristic argument showing the particular nature of the P — 1 fluctu- 
ations that deflnes the so called 1// noise. If S{f) oc 1/ and G{t) oc l/r'"', then the 
cosine Fourier transform above roughly shows that 1//^ has to be proportional to 1/ f^~'^. 
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Consequently, when (3 approaches one, 7 has to vanish. This means that for (3 = 1, the 
assumed form of G{t) is replaced by a logarithmic behavior implying a much lower rate 
of decreasing of the correlation i.e. correlation times grow strongly. 

Another case of particular interest is the exponent /3 = 5/3 corresponding to the 
homogeneous turbulence in fluids. The analogy is related to the presence of hierarchical 
processes like cascades from large to small scale fluctuations. 

The analysis of sunspot data along the lines just described provided an exponent 
clearly compatible with /5 = 1 as it is shown in figure 2. 
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Figure 2: Power spectrum of the sunspot number fluctuation. The arrow indicates the 
peak at / = (11 yr)-^ = 0.0076 month'^ 

In order to improve the analysis, a standard Butterworth filter [7] ( 1/[1 + cu^"] ) for 
n = 6 and 8, was applied to the data in order to separate the evident low frequency 
components of the signal under consideration from the high frequency ones. The low 
frequency part shows clearly the well known ~ 11 year period of activity as a sharp peak 
in the corresponding spectrum, while the high frequency component has the noisy behavior 
presented in figure 2. A power spectrum study of this high frequency part provides the 
result 

p{f) ^ r'-'^'-' (5) 

for 

1/(132 month) < f < (1/1 month) 

or alternatively 

1.3 X 10"^° Hz < f < 4 X 10"^ Hz (6) 
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We then conclude that the power spectrum of the sunspot number temporal fluctuations, 
along more than two hundred years, is well represented by a 1// behavior. 

3.2 Detrended Fluctuation Analysis (DFA) 

DFA [8] is an alternative statistical method to classify temporal correlations in a time 

series. It consists in dividing the time series of n values x{t) into n/s nonoverlapping 
groups, each with s values. Inside each group, the local trend or tendency y{t) = at + 6 is 
defined as usual through the linear least-square fit of the values in the group. From this, 
the detrended fluctuation function F{s) is defined as 

[ns)f = - E W)-y{t)f ; = 0, !,...,(- -1 (7) 

* t=ks+l ^* ^ 

Finally, the average of F{s) over the n/s intervals is computed. 

Whenever the time series data x(t) is randomly uncorrelated (or short range corre- 
lated), the expected behavior for the average is a power law, namely 

{F{s)) ~ (8) 

with an exponent Ha = 1/2, characteristic value of a random walk. If for a given range of 
s the exponent Ha happens to be different from 1/2, it detects the presence of long-range 
correlations in that range. When Ha > 1/2 there is a persistent behavior or a positive 
correlation between the increments, while if Ha < 1/2 the behavior is antipcrsistcnt 
showing the presence of a negative correlation. In both cases it is customary to speak, 
after Mandelbrot, about fractional Brownian movement. 

The DFA avoids inherent trends at all scales and allows to easily probe local correla- 
tions. Notice also that Ha is related to the fractal dimension through 

D = 2-Ha (9) 

if the distribution is described in a two dimensional space as is the case for the time series 
under analysis. 

The behavior of the detrended fluctuation function for sunspot data is shown in flgure 
3. We found the presence of a persistent behavior characterized by the exponent Ha — 
0.62. 

3.3 Jump-Size Distribution Method 

A further analysis of the time series can be performed in terms of the jump-size distribu- 
tion. The main objective is to determine the probability distribution of the logarithm of 
consecutive relative jumps in the number Si of sunspots, namely 

(10) 



Li = log 
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Figure 3: Average of the detrended fluctuation. 

together with the corresponding distribution 

Fl = ProbiL < L) (11) 

In figure 4, logroll ~ Fi{L)] is displayed for the case of the sunspot number using the 
data plotted in figure 1. It is found that Fl shows an inverse power behavior of the Zipf 
type [10], i.e., a straight line on a double logarithmic plot, showing that [1 — Fl{L)] can 
be expressed as the power 

1-Fl{L)=L-- (12) 

We have found two different regimes: The one corresponding to negative L with ti 0, 
and the one corresponding to positive L, characterized by 

r2 = 1.23 ± 0.20. (13) 

In the Appendix we have included an analysis of the probability distribution of relative in- 
crements, just to exphcitly show that this power behavior should be expected particularly 
in the case of a probability distribution of events following a power law. 

4 Final Comments 

It is well-known that the 1// noise is one of the most common and ubiquitous features 
in Nature. Superposition of effects giving rise to signals with scale invariant distributions 
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Figure 4: Logarithm of the cumulative distribution of consecutive relative jumps, logio[l — 
Fl{L)] (see equation (11)), plotted versus L for the case of the number of sunspots 

of correlation times, the so-called scale similarity, could be on the basis of the observed 
behavior. However, a proper explanation for such a behavior is still lacking and for that 
reason, the physical origin of the 1// noise is pretty much an open question. Nevertheless, 
a recent analysis [11] of models presenting self-organized criticality seems to indicate that 
the origin of the 1// behavior could be found in the superposition of local power-spectra 
with characteristic frequencies suitably distributed in space and having spectra that do 
not present the 1// behavior. It has been also shown, in the mentioned context, that 
one of the main features relevant for describing 1// noise is the well-known process of 
activation / deactivation. ^ 

We find that both mentioned ingredients are present in the case of sunspots. In fact, 
the noise in the magnetic field measured at earth, for short times (large frequencies) 
even during magnetic storms, presents a 1//^ spectrum. Consequently, the property of 
the sunspots spectrum can be related to the interference of random distributed magnetic 
disturbances. Moreover, the commonly accepted mechanism for the sunspot appearance 
is an activation/deactivation process, as it was discussed in the introduction. 

With reference to the roughness exponent, one can conclude that it also detects the 
persistent behavior or positive correlation between increments, due to the physical mech- 
anisms giving rise to sunspots. 

^For description of different examples of these avalanche type of processes, in connection with 1// 
noise, see for instance [11, 12, 13, 14]. 
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Appendix: On Probability Distribution of Relative In- 
crements 

In this Appendix we illustrate how the uniform and the power distributions give rise to 
the probability distribution of relative increments of the Zipf-like type. To start with, let 
us consider a random variable X with probability distribution fx{x) (—00 < x < 00). 
The distribution fx is normalized, that is. 



f 



fx{.x)dx = 1. (A.l) 
The cumulative probability distribution Fx{x) is defined via: 

Fx{x) = r fx{z)dz. (A.2) 



Fx{x) gives the probability that X < x, and clearly, fx{x) = dFx{x) / dx. Let Xj, i = 
1,2,.. ., be a series of independent instances of X, that is, each Xi is an independent sample 
of the distribution fx{x). We are interested in calculating the probability distribution of 
certain functions of these numbers, namely, rj, i?j and Lj, i — 1, 2, . . . defined as 

n = ^-1 (A.3) 

Xi+l 

Ri = \ri\ (A.4) 
Li = hiRi (A.5) 

The probabihty distribution of the r^'s can be evaluated as follows (overhned symbols 
refer to random variables): 

Fr{r) = Prob(r<r) (A.6) 
= Prob(^-l<r) (A.7) 

= Prob(a;i < (r + 1)0:^+1). (A.8) 

Since x.-. and independent, the probability of the last equation can be evaluated 

st r aight f orwar dly 

/oo 
fx{xi+i)Fx ((r + dxi+i. (A.9) 

-00 
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Or, equivalently, 

dxjxix,) / dx2fx(x2). (A.IO) 
The probability density distribution, can be then evaluated via: 

rrfa;a;/x(a;)/x((r + l)a;). (A.ll) 
ar J-oo 

The distributions for R and L can be calculated similarly: 

Fr{R) = Proh(R <R) = Prob(|r| < R). (A.12) 
It is evident that Fr{R) = for < 0. For it! > one has, 

Fr{R) = Prob(-i? <r<R)= Fr{R) - Fr{-R). (A.13) 

On the other hand, 

Fl{L) = Prob(I<L) (A. 14) 

= Prob(lni?<L) (A.15) 

= Prob(i? < e-^), (A.16) 

therefore, 

=F^(e^). (A.17) 

The probability density distribution can be then calculated similarly as in equation (A.ll). 

Let us now evaluate these probability distributions for the particular cases of the 
uniform and the power distributions. 

The uniform distribution 

The uniform distribution is given by 

r 1 if < X < 1 
fx{x) = (A. 18) 

[ otherwise 

Applying equation (A.IO), one obtains 

r if r < -1 



Fr{r) = < 



r + 1 

—— if - 1< r < 
2 - - 

2r + 1 .„ 

— if r > 

I 2(r+l) 



(A.19) 
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And then from equation (A. 13) 



r 



MR) = 



2R+1 1-R 



2{R + 1) 

2R+1 
[ 2{R + 1) 



if i? < 



if < < 1 



if > 1 



The power distribution 

Consider the distribution 



fx{x) = 



' (a - l)x-" if X > 1 
otherwise 



where a is a real parameter {a > 1). 

Applying equation (A. 10), one obtains 



r 



Fr{r) = < 



1 - 



1 1 



r + 1 2 



1 



2(r+ l)°-i 



And then from equation (A. 13) 



r 



Fr{R) = < 



if r < -1 

if - 1 < r < 

if r > 

if < 



-l)°-i ^ ^ V1--R 2j 



1 - 



2{R+1) 
1 

2(i?+ l)°-i 



if i? > 1 



(A.20) 



(A.21) 



(A.22) 



(A.23) 



In these two examples, two different straight line regimes are clearly present, as it is 
illustrated in figure 5 where 1 — Fl{L) (see equation (A. 17)) is plotted as a function of 
L. This is precisely the behaviour of the corresponding distribution in the case of the 
sunspot number, as shown in section 3.3 (figure 4). 
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Figure 5: Some cumulative functions (1 — Fl{L)) plotted versus L. (a) The uniform 
distribution (equation (A. 20)). (b) The power distribution (equation (A. 23)) in the case 
a = 3/2. (c) Same as (b) but for a = 3. In every case two asymptotic regimes are clearly 
seen. 
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